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Abstract 

The bifurcation diagram associated with the logistic equation u n+1 = at> n (l - u n ) 
is by now well known, as is its equivalence to solving the ordinary differential equation 
(ODE) it' = au(l — u) by the explicit Euler difference scheme. It has also been noted by 
Iserles that other popular difference schemes may not only exhibit period doubling and 
chaotic phenomena but also possess spurious fixed points. We investigate computationally 
and analytically Runge-Kutta schemes applied to both the equation u' = au(l - u) and 
the cubic equation u' = au(l-u)(6-u), contrasing their behaviour with the explicit Euler 
scheme. We note their spurious fixed points and periodic orbits. In particular we observe 
that these may appear below the linearised stability limit of the scheme, and, consequently 
computation may lead to erroneous results. 


f This work was performed whilst a visiting scientist at NASA Ames Research Center, 
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I. Introduction 

It is now well established that numerical schemes for solving ordinary differential equations 
(ODEs) exhi bit period-doubling and chaotic behaviour when used with time steps above 
their lineari sed stability limit. The most well know example of this is the explicit Euler 
difference scheme when applied to the ODE 

u' = au(l — u). (1-1) 

For this equation the scheme becomes 

u «+i _ u n^. aAtu n (l — u n ), (1-2) 

where At is the timestep being used. Figure 1 shows the bifurcation diagram obtained 
for the scheme. The bifurcation diagram is a plot of u n against r = aAt for two hundred 
iterates after first 600 iterates have been taken to allow the solution to settle. As can be 
seen, for r < 2, which is the linearised stability limit at the stationary point u = 1, all 
the successive iterates take the value 1, the stable equilibrium of the differential equation. 
Above this value of r the iterates alternate between two values whilst for even larger values 
of r the iterates cycle among four distinct values and so on. This phenomenon is known as 
period doubling, which for this case degenerates to chaotic behaviour where no finite set 
of distinct values can be discerned. Finally at r = 3 the numerical scheme breaks down 
as its solutions become attracted to the attractor at infinity. Notice however how the 
period-doubling behaviour is interupted by solutions of lower periods, a feature of most 
bifurcati on di agr ams of simple discrete maps [8]. The numbers labelling the branches of 
the bifurcation diagram indicate their penod (upto period 8), in addition the subscript E 
on the period one label indicates that it is an essential fixed point i.e. a fixed point of 
the ODE (1.1). We shall see later that other, spurious, fixed points may be produced by 
some numerical schemes. 

This type of period doubling behaviour is well known, the above example being equiv- 
alent, after a linear transformation, to the logistic equation of population dynamics [6] 

t> n+1 = ov"(l — v n ), (1.3) 

which is known for its chaotic behaviour. 

It is to be noted that once the explicit Euler scheme exceeded its linearised stability 
limit it announced the fact by its period 2 behaviour where the solution oscillates between 
two values. This is because linear multistep methods, of which the explicit Euler scheme 
is a simple example, have only the fixed points of the differential equation (Iserles [2]). 
This means that if the iterates take on a single value, i.e. a period 1 solution, then this 
is a solution of the differential equation. However Iserles [2] showed that for the class of 
Runge-Kutta schemes this need not be the case and, as we shall show later, these schemes 
can produce spurious solutions which are period 1 (non-osdllatory) but are not solutions 
to the differential equation. We investigate this phenomenom for three popular Runge- 
Kutta schemes and observe that such spurious solutions may, in some circumstances, be 
obtained for values of r below the linearised stability limit. We also note that for these 
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schemes we must sometimes greatly exceed the linearised stability limit before oscillatory 
behaviour hints towards a spurious solution which has been present since the limit was 
reached. Finally we observe that, even for the explicit Euler scheme, more than a single 
period-doubling solution may exist, the choice being dictated by the initial data, hence 
bifurcation diagrams become incomplete and fragmented as the iterates jump from one 
solution to another as r increases. We should note at this point that the explicit Euler 
scheme can, as well as being a linear multi-step scheme, be considered also as a first order, 
one-stage Runge-Kutta method. 

The implications of behaviour detailed above ranges far beyond pure ODEs. For many 
steady state calculations for partial differential equations (PDEs) numerical calculations 
are performed using ODE schemes, often Runge-Kutta, to ‘time march the solution. Could 
therefore the behaviour observed in this paper explain non-convergence experiences where 
the solution of the PDE appears to oscillate between more than one steady state, or where 
the residual will decrease only so far before reaching a plateau? Indeed, even though the 
mechanisms involved are far more complicated than those studied here, this could well be 
an exp lanat ion. This then leads us to ask if the solutions which are obtained, without any 
oscillatory behaviour to indicate to the contrary, are the true solutions to the differential 
equations? 

We do not try to answer these questions here, rather initiate studies which one day 
may allow us to supply such answers with confidence. We not only provide the numerical 
evidence of bifurcation diagrams, but trace the lower periodic fixed points of the schemes, 
where possible analytically. A description of the implication of the dynamical behaviour of 
finite difference methods for practical PDEs in computational fluid dynamics is described 
in Yee and Yee et al. [10,11]. Reference [11] can also be used as an introductory and state 
of the art of the current subject. A general account of the theory of asymptotic states of 
numerical methods applied to initial value problems may be found in Iserles et al. [3]. 

The schemes considered here are 

1) Explicit Euler scheme (for comparison) 

2) Modified Euler (2nd order Runge-Kutta) 

3) Improved Euler (2nd order Runge-Kutta) 

4) Heim’s scheme (3rd order Runge-Kutta) 

5) Runge-Kutta 4th order 

being tested both on the ODE (1.1) with quadratic forcing term and on the ODE 

u' = otu(l — u)(b — u), (1-4) 

for constant 0 < b < 1, with its cubic forcing term. 

In the next section we analyse the differential equations and the schemes, fin ding their 
fixed points, in Section III we investigate the local behaviour of bifurcations to spurious 
fixed points, whilst in Section IV we look at the higher order periodic orbits of the numerical 
schemes. 
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II. Fixed Points of the Equations and the Schemes 

Before we look in depth at the bifurcation diagrams for the schemes applied to the differen- 
tial equations we first look at some basic theory. We investigate the differential equations 
themselves, finding their fixed points, before introducing the s chem es, findin g their fixed 
points also. 

Consider first the general ordinary differential equation 

= /( u )> ( 2l1 ) 

where f(u) is a non-linear function of the variable u. 

This equation has fixed points u* (also know as equilibrium points, critical points or 
steady-state solutions) when 

/(«*) = 0 , ( 2 - 2 ) 

i.e. when the equation (2.1) is in equilibrium. If the fixed point is stable then u will be 
attracted towards it, otherwise if it is unstable then u will be repelled away from it. To 
discover the stablility, or otherwise, of a fixed point u* we must linearise the differential 
equation (2.1) about it. 

First we set 

u = u* + 5, (2.3) 

where we refer to 8 as a perturbation. Substituting this into the differential equation we 
obtain 


(«* + sy = f(u* + s), 

*/(«*) +W+-- ( 2 - 4 ) 

Since /(«*) = 0 and u m> = 0 we have the differential equation for 6 

8' = f'(u*)6 + 0(8 2 ), (2.5) 

which, ignoring the 0(8 2 ) term, has solution 

8 = tl'W. ( 2 . 6 ) 

We now see that if /'(u*) > 0 then the perturbation will grow and so the fixed point is 
unstable, whereas if /'(«*) < 0 the perturbation will not grow yielding a stable fixed point. 
If we have f'(u m ) = 0 we must consider a higher order perturbation (i.e. not ignore the 
higher order terms of 8 in (2.5) ) to determine the nature of the fixed point. 

If we turn now to the two differential equations we are considering, namely 

u' = mt(l — u) (2.7) 


and 


u' — mi(l — u)(b — tt), 


( 2 . 8 ) 


4 


we can find and investigate the stability of their fixed points. For (2.7) we have /(u) - 
au(l - u) which is zero at «* = 0 and u* = 1. The derivative of / is / (u) - oc{l- 2u) 
which takes the values +a and -a at these fixed points respectively. Therefore (2.7) has a 
stable equilibrium point at u* = 1 and an unstable equilibrium point at u* = 0, assuming 

that a is positive. _ 

For the second ODE (2.8) we have f(u) = au(l — u)(b—u) which has zeros at u — U,1 

and b. The derivative is 

f'(u ) = a[(l — u)(b — u) — u(6 — u) — u(l — u)] (2.9) 


which takes values aft, a(l - b ) and -ab(l - 6) at these fixed points respectively and so, 
since 0 < b < 1, we see that u* = b is a stable equilibrium point whilst u* = 0 and u* — 1 

are both unstable. . _ . 

We next look at the numerical schemes which we are investigating. First there is the 


explicit Euler scheme 


u n+l =u n + At/(u"), 


( 2 . 10 ) 


which is a linear multistep method. The remainder of the schemes are Runge-Kutta meth- 
ods, which are non-linear in f(u). The first of these is the second order Runge-Kutta 
method known as the modified Euler scheme and is given by 


u n+i = tt » + At/(u n +*A* /(«")) , 


( 2 . 11 ) 


whilst another second order Runge-Kutta method is the improved Euler scheme 

u «+i = „» + j Af [/(u n ) + / (u n + A t/(u n ))] . (2.12) 

The higher order schemes that we investigate are the third order Eeun scheme 

ti" +1 =u n + ^i+3A:3), 

*1 = /( O, (2-13) 

k 2 =/(« n +iA<A: 1 ), 

= f(u n +|A<fc 2 )t 


and the fourth order Runge-Kutta scheme 

u n+1 = u n + ^(jfcx + 2J fc 2 + 2kz + ^4), 
6 

*i = /(« n ), 

k 2 = f(u n +*At*i), 

^3 = /( u ” +3 Atfc 2 )i 

k*. = f{u n + A tk 3 ). 

If we consider a general scheme of the form 


(2-14) 


u n+1 =u n + AtF(u n ,At), 


(2.15) 



which encompasses all those mentioned above for a given f(u), then we can also find the 
fixed points of the scheme, i.e. u* such that F(u*,At) = 0. Here we use the term fixed 
point to m ean a fixed point of period 1 as opposed to a fixed point of higher period 
(periodic orbit). Note that now these fixed points depend on the additional parameter At. 
We may also investigate the stability of these fixed points in a similar manner to that used 
to investigate the stability of the fixed points of the differential equation. First we perturb 
the fixed point, writing u n = u* + £ n , and substitute into (2.15). After expansion of the 
resulting term on the right-hand side, cancelling of the u* and neglect of high order terms 
in 6 n we obtain the difference equation 

S n+1 = s n [l + A *F«(u*, A*)] (2.16) 

for S n . This has solution 

6 n = [l + AtF u (xi m ,At)] n 6° (2.17) 

and so for stability we require 

|l + AtF tt (u*,At)| < 1, (2-18) 


—2 < AtF u (xi*,At) < 0. (2.19) 

Notice again that the parameter At appears (maybe not in a linear fashion) and so we 
shall speak of stability regions for certain ranges of At. 

Looking at the specific case of the explicit Euler scheme applied to the first ODE (2.7) 

we have 

u n+l _ u n + oAtu n (l — u n ), (2.20) 

i.e. F(u,At) = f(u) = au(l - u). This therefore has fixed points u* = 0 and 1, the same 
as for the differential equation. As noted by Iserles [2] this will always be the case for 
linear multistep methods. If we now investigate the stability of these fixed points we have 
F u = a(l - 2u) and so (2.19) gives that, for At > 0, u* = 0 is unstable, whilst u* = 1 
is stable for 0 < aAt < 2. This is known as the linearised stability range of the scheme 
applied to the differential equation, and we must choose a At within this range if we are 
to compute the correct stable solution. For the ODEs and schemes we are considering the 
ci may always be combined with the timestep At, acting as a scaling, and so for notational 
brevity we shall from now on use 

r = aAt. (2-21) 

The above stability region then becomes 0 < r < 2. 

If we repeat the process with the second ODE (2.8) and the explicit Euler scheme 
we obtain 0, 1 and b as the fixed points, the first two of which are unstable for r > 0 
whilst the latter (the correct stable equilibrium of the differential equation) is stable for 
0 < r < So, for example in the symetric case b = £ we must have 0 < r < 8. 

We may repeat this process for all the combination of schemes and equations, taking 
welcome aid from an algebraic manipulation package such as MAPLE [5] or DERIVE [1]. 
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The stable fixed points for the various schemes applied to equation (2.7), together 
with their stability ranges are given in Table 2.1, assuming that r > 0. The question 
marks indicate where stable fixed points are known to exist from numerical experiments 
but no closed analytic form has yet been found. 


Scheme 

fixed points 

stability 

Explicit Euler 

1 

0 < r < 2 

Modified Euler 

1 

0 < r < 2 


1 + 7 

0 < r < — 1 + y/5 w 1.236 


2 

r 

2<r<l + \/5w 3.236 

Improved Euler 

1 

0 < r < 2 


2-f*rdb\/r 2 — 4 
2r 

2 < r < >/8 « 2.828 

Heim 

1 

0 < r < 1 + (\/l7 + 4) » - (Vrf + 4)-^ « 2.513 


? 

? 

R-K 4 

1 

0 <r<|+(^ + fv'29)i 



+ (M_i v ^9)J a, 2.785 


? 

7 


? 

7 

Table 2.1 Numerical Period 1 Fixed Points of ti' = u(l — u) 


Scheme fixed points 

Explicit Euler \ 

Modified Euler 2 

l±-v/l-8/r 

2 

Improved Euler § 

i±\/ 1— 8 / r 
2 

1 v' 1 ~ 12 / r + v / ' 1 + 4 / r 

2 4^4 

1 , ■ Vl+4/r 

2 ’ 4 - C 4 

Heim i 


R-K 4 


1 

2 
? 

? 


stability 
0 < r < 8 
0 < r < 8 

8 < r < 4(1 + y/Z) « 10.928 
0 < r < 8 

8 < r < 12 

12 < r < 4(1 + V6) « 13.798 

12 < r < 4(1 + \/6) « 13.798 
0<r <4(l + (\/l7 + 4)i 

— (y/Vi + 4)~» ) « 10.051 

0 < r <? 

? 

? 


Table 2.2 Numerical Period 1 Fixed Points of u' = ti(l - u)(| - u) 
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The stable fixed points for the various schemes applied to equation ( 2 . 8 ) for the special 
symmetric case 6 = 5 , together with their stability ranges are given in Table 2.2, assuming 
that r > 0. The more general case is illustrated by the various Figures (see below). Again 
question marks indicate where the algebra has defeated us and the algebraic manipulators 
at our disposal due to the high order polynomials in two variables involved. 

Figures 2-4 illustrate these fixed points, as well as higher order periodic orbits of the 
schemes. Again the numeric labelling of the branches denote their period, although some 
for period 4 and 8 are omitted due to the size of the figures. The subscript E on the 
period one branch indicates the essential fixed point of the differential equation whilst the 
subscript S indicates the spurious fixed points introduced by the numerical scheme. 

The period one orbits in the figures where obtained by solving F(u*, At) = 0 numer- 
ically and checking a discretised verison of (2.18), the higher period orbits were obtained 
by numerical approximations of the conditions of Section III. 

As can be seen, except for the explicit Euler scheme (as expected) and the Heun 
scheme, there are spurious stable fixed points as well as the correct stable fixed point 
(ti* = 1 for (2.7) and u* = | for the symmetric case of (2.8)). Although in the majority 
of cases these occur for values of r above the linearised stability limit this is not always 
the case. In particular for the modified Euler scheme applied to (2.7) we see that there is 
stable spurious orbit below the linearised stability limit of r = 2 for u* = 1. This is outside 
the interval 0 < u < 1 and so it is unlikely that it will be picked up accidentally since 
usually initial value tt° would be chosen between the two fixed points of the differential 
equation. The fourth order Runge-Kutta scheme applied to the same equation, however, 
exhibits a spurious critical path which not only lies below the linearised stability limi t but 
also in the region between the fixed points of the differential equation and so could be 
easily achieved in practice. For more complicated problems such spurious points could be 
computed and mistaken for the correct equilibrium. 

Another dy namical behaviour of these spurious fixed points generated by the schemes 
is that when all such spurious paths are plotted for a given combination of equation and 
scheme they often resemble period doubling bifurcations.On the main branch, where period 
\e lies, the spurious paths are branching from the correct fixed points as they reach the 
linearised stability limit, and sometimes even forking again as r increases still further. 
The result of this is that bifurcation diagrams calculated from a single initial condition tx° 
will appear to have missing sections of higher period orbits (see Figures 8 - 11 ), or even 
seem to jump between br anch es. This is in fact not the case since all attractors of higher 
period orbits must be present, although as we shall see such higher period orbits may be 
non-unique, even for the explicit Euler scheme, thus propagating this effect throughout 
the bifurcation diagram. In order to compute ‘full’ bifurcation diagrams we must overplot 
a number of diagrams obtained using different starting values u°. For the higher order 
schemes many such overplots will be needed to fall within all the basins of attraction. 
Such diagrams are illustrated in Figures 5-7, their earlier stages resembling the fixed 
point of Figures 2-4. The terms tram critical and supercritical bifurcations [9] 

refer to the nature of the bifurcation of the fixed point as it reaches its linearised stability 
limit. Supercritical bifurcations have both of their branches being stable at the bifurcation, 
whilst for transcritical bifurcations one branch is stable whilst the other (at least initially) 
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is unstable. As can be seen those solutions with spurious fixed points below the linearised 
stability limit of the scheme are a result of transcritical bifurcations. 

In the next section we make use of perturbation arguments to investigate the local 
nature of the bifurcations from an essential fixed point of the differential equation to 
spurious period 1 rest states of explicit Runge-Kutta methods of order < 4. 
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TTT Bifurcation to Spurious Period One Solutions 

In th i g section we investigate the behaviour of explicit s-stage Runge-Kutta methods 
of maximum order in the neighbourhood of a stable fixed point u* of (2.1)* It is convenient 
to express the general s-stage method in the form 

u n+1 = u n + At ^ Cjf(zj), (3.1) 

j=i 

where 

Zj = u" + At 53 bj'ifizt), j = 1, 2, . . . , s (3.2) 

l=i 

and, in the more standard notation of the previous section, kj = f(zj). In order that we 
may assume the method to have order s, we s h al l have to restrict attention to the case 
s < 4 [ 4 ], One of the conditions that is necessary for the method to be consistent is that 

c T e = 1 (3.3) 

where c = (ci, c 2 , . . . , c,) T and e = (1, 1 , . . . , 1) T . For a general discussion of bifurcations 
of maps of the interval we refer to Whitley (9). 

Clearly u n = tz* is a fixed point of (3.1), (3.2) if u* is such that /(u*) = 0. Expressing 
the ma pping (3.1), ( 3 . 2 ) in the form (2.15) and linearising about u* we find that 

AtF„(u*, At) = ps? (I - pBy'e. (3.4) 

where p = Af/'(u*) and B denotes the s x s array of weights bjj that occur in (3.2). With 
$ < order s may be achieved with s stages and it is a well established result that the 
Jacobian of the mapping is given by 

1 + AtF u (u*, At) = 1+ p + + • • • + ^p\ 

an 0(p-») approximation to e p . Hence, 

£ T (I - pB)~ 1 §.= 1 + 1* (3*5) 

the first s terms in the Taylor expansion of (e p - 1 )/p. 

When u* is a stable fixed point of (2.1) we have /'(«*) < 0 and the linearised stability 
condition (2.19) will be satisfied for all At sufficiently small since, from (3.4) and (3.3) 

AtFUu*, At) = A t/'(u*) + 0(At 2 ). 

Hence, u* will be a stable fixed point of (3.1), (3.2) in some interval At € (0, At*] or, 
equivalently, p G [/>*, 0) ( p* = p*(s) = At* /'(«*) < 0). The interval [p*,0) is usually 
called the interval of absolute stability of the Runge-Kutta method defined by (3.1), (3.2). 
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It is easily shown that p* satisfies 


p'c T (I - p*B)~ 1 e = -2 


for s = 1 and a = 3 so that, by (3.3) and (2.16) 6 n+1 « -6 n as p decreases beyond p* 
leading to a period doubling (flip) bifurcation at p = p*. This situation will be described 

briefly in the next section. _ _ 

In the cases s = 2 and s = 4 it may be shown that the limit of absolute stability, p , 

satisfies . . 

p'c T (I - p m B)~ l e = 0 (3-6) 


so that |£ n+1 1 > |<$ n | for p < p* and there is a loss of stability of the fixed point u* that 
would, for a linear problem, lead to |u n | — ► oo as n — * oo. One of our aims here is to 
show that such divergence does not occur for a genuinely nonlinear differential equation 
constant) but that there is a bifurcation to a fixed point, u, that is spurious in 
the sense that f(u) ^ 0 so that it is not a stationary point of (2.1). We also indicate how 
the nature of the bifurcation is influenced by properties of the method and those of /(u). 

Any fixed point, u, of (3.1), (3.2) must satisfy 


5^c i /(z i ) = ° 
>=i 


(3.7) 


and 


Defining 


i-i 


Zj = xi + A t 6j,//(zj), j — 1, 2, ...» s. 


(3.8) 




e = tt — it*, (3*9) 

we seek a solution of (3.7) and (3.8) for p close to p* (At close to At ) in the form 

At = At* +ae-\-be 2 H (3.10) 


and 

Zj =U* + OLjt + 0j€ 2 + 7jC 3 + • • ' 

Prom this we deduce that 


j — 1, 2, . • . , s. 


(3.11) 


f(z s ) = </'(«>> + y [/»? + 2 /*(«•)&] 


+ € J 


!/'"(«*)«? + /"<«>>* + /vbi 


(3.12) 


+ 


Substituting (3.10)-(3.12) into (3.8) and equating like powers of e leads to 

(I — P*B)ql = e, 


(3.13) 
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and 


(3.14) 


(I - p'B)l = B 


a/'(«’)a+5Afr(“*)a a 


where a is to be determined and o = (c*i, • * • » and denotes the vector whose com- 
ponents are the squares of the corresponding elements of a. It is convenient for subsequent 
manipulations to rearrange (3.14) to read 




«fV)a+ . 


(3.15) 


Since (3.6) and (3.13) imply that sF'si. = 0 we find, on combining (3.7) and (3.12) and 
neglecting terms in e 3 , that the condition for a fixed point becomes 

/'V)cV+ 2/V)c t £ = 0 


or, with p * = f'(u*)At *, 

A f/'VWV + 2p*z T @_ = 0. 
Talcing this together with (3.15) we obtain 

At*/ M (u*) £ T (J — p*B)~ l o? 
a ~ 2f'(u m ) sl t (I - p m B)~ 1 QL 


(3.16) 


where a is given by the solution of the system (3.13). Thus a is well-defined provided that 
c T (I — p*B)~ 1 o. does not vanish. 

When a/0, that is, when 

/"(«*) ^ 0 (3.17) 

and 

£ r (J-p*B)-V^0 (3.18) 

there is a tr ans critical bifurcation given, to first order, by (3.9) and (3.10). The first of 
these conditions depends on the differential equation (and is violated, for example, for 
/(«) = au(l - tt)(l /2 - it) with u* = 1/2) and the second condition depends solely on the 
Runge-Kutta method. 

To study the stability of the bifurcating solution we write 


A(e) = 1 + A tF„(u, At) 


to denote the Jacobian of the mapping at u = u* + e, At = At* + ae H . We then find 

that . 

A'(0) = ^At , /"(u*)£ T (I — p*B)~ 1 s? 

and, since A(0) = 1, we shall have A(e) < 1 provided 

€/"(u*)c r (J - p*B)~W < 0. (3.19) 
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to a stable branch. Using (3.13), it 


This condition determines the sign of e that gives rise 
may be shown that 


c^(I — p*B) 1 a = — (1 + pcF(I — p m B ) e)|p=p* 

dp 

and, since the argument on the right is a decreasing function of p by the definition of 
interval of absolute stability, it follows that c T (I - p*B)~ 1 a < 0. This result together with 
(3.16) imply that the stable branch emanating from (tx*, A<*) is given by (3.9) and (3.10) 
for at > 0, that is A t > At*. 

If we consider the most general second order two stage Runge-Kutta method, it may 
be parameterised by 0) so that c = (1 — 6, 9')'^ and 


B = 


( 


0 0 \ 
1 /( 2 *) 0 ) 


It is easily deduced that p* = —2, qF ( I — p*B ) 1 ql — —1 for all 6 and 

SL T (I ~ P*B)~ 1 e^ = 

Thus, the method will not generate a transcritical bifurcation if 0 — 1/2 (the improved 
Euler method (2.12)). For 6 1/2 we obtain 

A **/"(«*) 1-20 
a ” /'(u*) “ 0 ■ 

In partic ular , for the modified Euler method ( 6 — 1) with f (u) = o;tt(l tx) and tx 1 
we obtain a - -At*. Thus, from (2.21), (3.9) and (3.10), the bifurcation occurs at r = 2 
and is described to first order by 

r « 2(2 — tx) 

and, for stability, r > 2 so that u < 1. These results are seen to agree with the graphical 
results shown in Figure 2. 

For the fourth order Runge-Kutta method, the limit of absolute stability is p* ~ 
—2.785 as given by the negative retd root of (3.5) with s = 4. Following a tedious calculation 
we find that (3.16) gives 

a « — 0.9598/" (tx*)/[/'(tx*)] 2 

so that a transcritical bifurcation always occurs provided /"(tx*) ^ 0. The equation of the 
tang ent line at the bifurcation point is given by 

tx « 1 + 0.521(r - 2.785) 

for /(tx) = mx(l - tx) with r > 2.785 so that tx > 1 for stability. These findings agree with 
the graphical results shown in Figure 2. 
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In those cases where either (3.17) or (3.18) is not satisfied we find that a = 0 in (3.10) 
and higher order expansions are necessary to determine the nature of the bifurcation. 
Omitting the details, we determine the coefficient of e 2 in (3.10) to be 

b= El L ~ elEEll**. (3.20) 

6/'(u*) 3 c T (I-p*B)~ 1 a 


where 


a* = diag(a)[(/'(ti*)/"V) - 3/"(u*) 2 )& 2 + 3/"(ti*) 2 (I - p*B)~ 1 o^]. 


When a = 0 and 6 # 0 the bifurcation at At = At* will be of limit point or pitchfork 
type (see Whitley [9]). Moreover, it will be supercritical if b > 0 and subcritical if fe < 0. 
Generally speaking the former will be stable and the latter unstable. 

Of the examples we have looked at in this section, only the improved Euler scheme 
leads to a pitchfork bifurcation through violation of condition (3.18). In this case we find 


that 


f(u*)/ w (u*)- 3/»(u«) 2 
3/'(u*) 3 


(3.21) 


The other instances of pitchfork bifurcations occur through failure of (3.17) and the ex- 
pression for b simplifies to 


6/'(u’) 2 £ r ( I — p*B) _1 SL 


This leads to /»'(„•) 1-3* + 30* 

b ~ /'( u *) 2 30 s 

for the general second order Runge-Kutta method and 


b = 2.187 


no 


/'( u *) 2 


( 3 . 22 ) 


for RK4. The factor involving 0 in (3.22) is always positive and it is interesting to note 
that the value obtained for b is the same for both 0 = 1 (modified Euler) and 0 = § 
(improved Euler). This is in accordance with the results shown in Figure 3 (at At = 9, for 
example). Thus, a (stable) supercritical pitchfork bifurcation results in all cases provided 
that f ,n (u*) > 0. For instance, the function f(u) = tz(l-u)(l/2-u) = -i(w-|)+(u-|) 3 
has three real zeros and, with u* = |, /'"(«*) = 6. On the other hand, for the function 
f(u) = — j(ti — 1) — (u — |) 3 which has only one real zero, u* = 5, f m (u*) = —6 and an 
(unstable) subcritical pitchfork bifurcation would result. 

The perturbation analysis described in this section has shown that it is possible to 
predict not only the onset of instability at an essential stationary point u* of the differential 
system but also to determine the nature of bifurcation that occurs and the stability along 
the bifurcating branch. In the next section we investigate some of the higher order orbits of 
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the schemes (a feature not present in the original differential equations which we consider) 
where consecutive iterates of (2.15) oscillate between two or more values. It was shown 
earlier in this section that such bifurcations occur for odd order methods. The orbits in 
these cases are generally much more difficult to obtain analytically, even with the aid of 
algebraic manipulation software. However, such analysis is presented where possible and 
numerical backup used for the harder cases. 
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IV. Periodic Orbits (Fixed Points of Higher Period) 

If we use a difference scheme to solve an ODE using a value of r which is slightly above 
the stability limit for a fixed point of the scheme (either spurious or one belonging to the 
differential equation) then we often find that the iteration (2.15) will oscillate between two 
values. This is known as period doubling, a process which is often repeated again at the 
stability limit for the period two orbit and so on. It is this process which is illustrated 
by the bifurcation diagrams. As well as period doubling, embedded regions of lower, even 
odd, periods often occur, frequently as a prelude to chaos. Chaos is the state where there 
is no finite set of attractors which the iteration visits. Finally instability of the scheme 
will usually set in when the iterates are attracted to the global attractor at infinity. These 
latter features are beyond the scope of this paper. We do however investigate the low 
order periodic orbits of the iteration (2.15) and illustrate that, like the fixed points of the 
scheme, non-unique stable periodic orbits may co-exist for given values of r. The orbit 
achieved in these instancies will depend on the initial value of the iteration u°. In such 
cases bifurcation diagrams produced using a single tx° will appear to have parts of their 
higher order orbits missing, whereas the true situation is the selection of only one of the 
possible lower period orbits. Bifurcation diagrams therefore can be data dependant, the 
full dingr*™ only being obtained by the super-position of several such ‘sub’ diagrams. 

Consider the case of period two orbits. This means that there must exist two values 

tz° and ti* such that _ . . . . 

ti = u + A tF{u , At) (41) 

ti* = u° + AtjF(ti°, At), 


that is, 

F(u°, At) + F(u*, At) = 0. (4.2) 

Rewriting this in terms of just one of the values, ti* say, we have that period two orbit 
states are given by the solutions of the equation 

F(u*, At) + F(u* -I- AtF(u*, At), At) = 0. (4.3) 


There are likel y to be more than two solutions of (4.3) and so they must be paired using 
(4.2). Note that the fixed points of the scheme will also satisfy both (4.2) and (4.3). 

To investigate the stability of the period two orbits we again perturb them slightly. 
For schemes of the form (2.15) which we are considering this amounts to the following (see 
e.g. Sleeman et al [7] for a technique applied to linear multistep schemes). Perturb the 
orbit state to u rt_1 = it* + <5 n_1 then linearising 

u »+i _ u »-i + , At) + F(u n_1 + AiF(u w “ 1 , At), At)] (4.4) 


about u* gives 

6 n+i = gu-ifr + AtF tt (u*,At)][l + AtF u (u* + AtF(u*, At), At)] (4.5) 

and so we see that for stability of the orbit we require 

|[1 -l- AtF u {u\ At)][l + AtF„(u* -I- AtF(ti*, At), At)]| < 1. (4.6) 
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These calculations are too involved for the higher order schemes that we are consider- 
ing, however we present the analytic forms for the period two orbits of the explicit Euler 

scheme in Table 4.1. 


Equation 
u' — att(l — u) 

u' = au(l — u)(j — ti) 


period 2 orbits 

r+2±'/^4 

2r 

ldVl-8/r 

2 

1 y/l—12/r x ■y /l+4/r 

2 4 X 4 

1 y/l-njr . y/l+Zjr 

2 + 4 X 4 


stability 

2 < r < -\/6 w 2.4495 
8 < r < 12 
12 < r < 14 
12 < r < 14 


Table 4.1 Period 2 Orbits of the Explicit Euler Scheme 

We note two things. Firstly that the period two orbits for the explicit Euler scheme 
are exactly the spurious fixed points of the improved Euler method, although their stability 
range is different. This is not too surprising when one considers that (4.3) with F(u, At) = 
/(tx), as is the case for the explicit Euler scheme, is precisely the equation for the fixed 
points of the improved Euler scheme. Secondly we note that for the second equation, (2.8), 
there are three different period two orbits, two of which co-exist over a range of r. We 
remark on this to point out that although linear multistep schemes possess only the fixed 
points of the equations their higher period orbits can be non-unique. 

Although we do not present analytic forms for the orbits (except those above), the 
orbits, upto period 8, depicted in Figures 2-4 were obtained by numerically solving (4.3), 
or the appropriate higher order form, and checking a discretised form of the appropriate 
stability condition, for example (4.6). 

Finally complete bifurcation diagrams for the various combinations of s chem es and 
equations, including cases of (2.8) where 6 ^ are shown in Figures 5-7, whilst bi- 
furcation diagrams produced using only a single initial data tt° are given in Figures 8 - 
11 illustrating the apparent missing of branches described above, together with the corre- 
sponding full bifurcation diagrams obtained by overlaying multiple single data diagrams. 
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V. Summary 

We have investigated the fixed points and periodic orbits of four Runge-Kutta schemes, 
contrasting them with those of the explicit Euler scheme which is known to possess only 
the fixed points of the differential equation. We have seen how not only do these schemes 
produce spurious fixed points but that these spurious features of the schemes can manifest 
themselves below the linearised stability limit for the correct fixed points. This raises the 
possibility of erroneous results when such schemes are used for computations on problems 
where the correct result is not known a priori. We have also observed how multiple orbits 
of a given period may co-exist, the particular one selected by the scheme depending on the 
initial data. Thus bifurcation diagrams produced using a single starting value may appear 
to be mis sing branches of higher order orbits. 

Future work will be directed towards investigation into the effect of using such ODE 
solvers for the source term component of reaction-convection equations. 
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Figure 2. Stable fixed points of period 1,2, 4,8 for u' — aru(l u). 
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Figure 5. Bifurcation diagrams for u' = atx(l — u). 
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Figure 9. Bifurcation diagrams (transcritical) for Modified Euler (R-K 2) scheme applied to 
u' = atx(l — u)(0.4 — u). 
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